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Abstract 

We provide a detailed analysis of 12 XMM observations of the nearest persistent 
extragalactic ultraluminous X-ray source (ULX), M33 X-8. No significant spectral 
evolution is detected between the observations, therefore we combine the individual 
observations to increase the signal-to-noise ratio for spectral fitting. The combined 
spectra are best fitted by the self-consistent p-free disk plus power-law component 
model with p = 0.571+°ofsl, kT,,, = l.SStofs keV, and the flux ratio of the p-free disk 
component to the power-law component being 0.63:0.37 in the 0.3 - 10 keV band. 
The fitting indicates that the black hole in M33 X-8 is of ~ IOM0 and accretes at a 
super-Eddington rate (~ 1.5 LEdd), and the phase of the accretion disk is close to a 
slim disk [p = 0.5). We report, for the first time, that an extra power-law component 
is required in addition to the p-free disk model for ULXs. In super-Eddington cases, 
the power-law component may possibly result from the optically thin inner region of 
the disk or a comptonized corona similar to that of a standard thin disk. 
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1. Introduction 

Ultraluminous X-ray sources (ULXs) are point-like, non-nuclear X-ray sources with 
isotropic luminosities of about 10^^ — 10^^ ergs s ~^ in nearby galaxies (Fabbiano 1989). These 
objects are interesting since their luminosities are intermediate between the luminosities of 
Seyfert galaxies (Lx ~ lO'^^ — 10^^ ergs s ~^) and those of black hole X-ray binaries (BH XRBs; 
typically Lx ^ 10^^ ergs s ~^). If X-rays are emitted isotropically below the Eddington limit, 
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the intermediate luminosity would indicate that ULXs harbor intermediate-mass black holes 
(IMBHs; Mbh ~ 20-10^ Mq). However, IMBHs are not required if the emission is anisotropic 
or relativistically beamed (King et al. 2001; King 2009; Kording, et al. 2002), or the accretion 
is super-Eddingtonian by a factor of a few (Begelman 2002; Poutanen et al. 2007). 

X-ray spectral fitting is one of the best methods known to weigh the black hole in ULXs 
while the dynamical method is unavailable. The sum of a multicolor disk (MCD; Shakura & 
Sunyaev 1973) and a power-law (PL) model is widely used to describe the X-ray spectra of black 
hole binaries. This canonical MCD+PL model reflects the expectation of the thermal emission 
from a standard thin disk around a Schwarzschild black hole along with the hard emission from 
the inverse Compton scattering of disk photons. Fitting ULX spectra with the same canonical 
model often shows that the PL component dominates the 0.3-10.0 keV spectrum and the disk 
component characterizes the feature of the soft spectrum significantly below 1 keV (Miller et 
al. 2003; Gongalves & Soria 2006). The disk temperature is inversely related to the black hole 
mass (Teff oc M-i/^ e.g., Makishima et al. 2000) in a MCD model, i.e., a heavier black hole 
tends to accrete with a cooler disk. The obtained disk temperature around 0.15 - 0.2 keV 
in many ULX spectra was used to suggest the existence of IMBHs with masses of ~ lO^M© 
(Kaaret et al. 2003; Lorenzin & Zampieri 2009). 

Meanwhile, hotter disks with temperature /cTin ~ 1-2.5 keV were also reported in lit- 
erature. In these cases, the disk component dominates the X-ray emission (Stobbart et al. 
2006), in a way similar to Galactic BH candidates. The hot MCD model (HD model) with fit- 
ted temperature fcTin ~ 1-2.5 keV would imply super-Eddington luminosities and stellar-mass 
black holes in the ULXs. However, the MCD model is based on the standard thin disk model 
and is valid only for luminosities well below LEdd (^ 0.1 ^Edd)- It is known that the radial 
temperature of a standard thin disk follows a power-law form as kT oc r~° ''^. Such a form is 
derived from the assumption of energy balance between the viscous heating and the radiative 
cooling. For sufficiently high accretion rates, the radiative cooling itself can not balance the 
viscous heating and the advective cooling becomes important or even dominant. In this case, 
the power-law form kT oc r~° ''^ is invalid. In fact, there are observed deviations from the stan- 
dard disk spectrum in BH XRBs as they approach their Eddington limit (Kubota & Makishima 
2004). In the slim disk model, that is applicable to super-Eddington accretion (Abramowicz et 
al. 1988), the radial temperature follows kT oc r~°'^. This demonstrates that the HD model is 
not self-consistent for ULXs (Gongalves & Soria 2006). 

The so-called p-fiee model also obeys a power-law form as Tcs oc r~^, but the temperature 
gradient p is allowed to vary from 0.5 to 0.75 (Mineshige et al. 1994). The p-free model is an 
extended MCD model and is a combination of two types of optically thick disk models, namely 
the standard thin disk and the slim disk. Vierdayanti et al. (2006) successfully applied the p-free 
model to the spectra of four ULXs, which were previously reported to contain cool disks and 
thus IMBHs. The fitting yields p values of ~ 0.5, consistent with the slim disk model, suggesting 
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that the black holes in these four ULXs have stellar masses and accrete at super-Eddington 
rates instead. Note that in their fitting, no extra PL component was considered. 

The ULX source M33 X-8 with an X-ray luminosity of ~ 10^^ ergs s~^ was discovered by 
Long et al. (1981) with the Einstein satellite. Although the position of M33 X-8 coincides with 
the optical center of the galaxy (La Parola et al. 2003), the hypothesis of an active galactic 
nucleus (AGN) is inconsistent with the estimated upper limit of 1500 M© on the central black 
hole mass in M33 (Gebhardt et al. 2001); moreover, no AGN activity has been found in other 
bands for this source. The source is of particular interest for many reasons. First, its X-ray 
spectra apparently prefer the HD model, but the model is unlikely to be appropriate at this 
luminosity. Second, it is the nearest persistent extragalactic ULX (Foschini et al. 2004). Third, 
it had up to 12 observation data available in the XMM- Newton Public Archive. In this paper, 
we present a detailed analysis of the 12 XMM exposures on M33 X-8, especially the combined 
spectra, which have much higher signal-to-noise (S/N) ratio and enables us to provide further 
constraints on the nature of this source. 

2. Data Reduction 

In Table 1 we list all the 12 XMM observations of M33 X-8, which were obtained from 
August 2000 to July 2003. Hereafter we refer them as Obs #1 through Obs ^^12 for convenience. 
The data were reduced with the XMM-SAS software version 7.1.0. To exclude intervals with 
background flares, we created light curves for photons above 10 keV, and used a count rate 
cut-off criterion to filter the light curves. The exact value of the cut-off was allowed to vary from 
field-to-field, to provide the best compromise in each case between excluding high background 
periods and facihtating the longest available exposure (Stobbart et al. 2006). We selected the 
data from good time intervals, by setting FLAG = and PATTERN < 4 for PN data, and 
PATTERN < 12 for MOS data. The source spectra were extracted from circles with radius of 
35" and centered at the nominal position of M33 X-8 (RA = 0l''33""50f89, Dec = +30°39'37!'2, 
J2000), while the background spectra were extracted from the same CCD chips as the source 
and at a similar distance from the readout node. The high spatial resolution Chandra image 
has confirmed that there is no obvious contamination from neighbouring point sources within 
several arcminutes (Dubus et al. 2004). For Obs #1 that the MOS camera was operated in 
small window mode, we used the background in the closest chip. With the SAS task epatplot, 
we found that three observations (Obs ^^2, #4, and #8) were affected by pile-up. Only PN 
data for these three observations were used and the spectra were extracted in annulus regions 
with radius 15" and 40" to circumvent the effects of pile- up (Loiseau 2005). In these cases, 
ARE files were calculated to correct the missing part of PSF, thus the correct flux level could 
still be measured through spectra fitting. In several cases, PN or MOS data were unavailable 
because the source was not covered by the detector, or due to CCD gaps. When available, 
PN and MOS data from each individual observation were fitted together for spectral analyses. 
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The 0.3 - 10.0 keV band spectra were fitted with the HEAsoft X-ray spectral fitting package 
XSPEC 12.3.1. All spectra were rebinned to have at least 20 counts per bin to enable the use 
of statistics. 

3. Spectral fitting 

3.1. Individual observations 

We first fit the 12 individual spectra with an absorbed MCD+PL model (diskbb+po 
in XSPEC) and a p-free disk model (diskpbb), respectively. We get similarly good fits for 
both models (Table 1). During the fitting, the absorption column density is allowed to vary 
(Feng & Kaaret 2006). Due to the limited number of photons, it is difficult to judge whether 
the MCD+PL model or the slim disk model is preferred. However, by examining the fitting 
residuals of the p-free model, in three of the observations with the highest number of data bins 
(Obs #1, #5, and #12), we find a weak hard tail above 7 keV, which we will further investigate 
for the composite spectra in §3.2. 

M33 X-8 was reported as a persistent source by Foschini et al. (2004), and a small flux 
variation was detected previously for a modulation of ~ 20% with a period of 106 days (Dubus 
et al. 1997). Consistently, we obtain that the amplitude of the flux variation between the 
observations is at < 20% level, except for Obs #4, and #8, which are about 1.5 - 1.7 times 
brighter than the average level of the other 10 observations. In Figure 1 we plot the X-ray light 
curve for the 12 observations of M33 X-8. The 0.3 - 10.0 keV luminosities are taken from the 
best-fitted MCD+PL model. While different models yield slightly different luminosities, the 
general pattern of the light curve will not be changed. The mean luminosity (with Obs #4 and 
#8 excluded) along with the ratio of the luminosities to the mean value are also shown. From 
the figure it is seen that M33 X-8 remains persistent during most of the observations, and we 
detected no significant spectral evolution through spectral fitting (see Table 1). No significant 
rapid variation was detected either within individual observations. 

3.2. Combined spectra 

Inspired by the fact of no significant spectral evolution between the individual observa- 
tions, we combined all the observations except for Obs ^^4 and #8, and obtained one co-added 
PN spectrum, one co-added MOSl and one M0S2 spectrum, respectively, in order to increase 
the S/N ratio of the spectra. The co-adding was performed with the FTOOLS "addspec" 
which adds pulse-height amplitude (PHA) spectra and background PHA files; detector redis- 
tribution and ancillary response were also combined with source net photon counts detected 
in each observation as co-adding weight. We fit the co-added PN, MOSl and M0S2 spectra 
simultaneously. We excluded Obs #4 and #8 because they are at significantly higher flux level, 
however, we have checked that including them does not change our results in this paper. 

We apply the same models mentioned in §3.1 to the combined high S/N spectra to test 
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whether these models can still provide adequate fits. Neither the HD model nor the simple 
p-free model can give good fits below ~ 1 keV (see Figs. 2a and 2b). La Parola et al. (2003) 
found that a thermal plasma component is required to represent the extended emission around 
the point source (also see Schulman & Bregman 1995). Accordingly, we add a Raymond-Smith 
component to the above models, then the is reduced by more than 50 with 3 additional 
free parameters, yielding an F-test probability lower than 10^^. After including the Raymond- 
Smith component, both models (HD and p-free) yield similar x^, and the fitting results are also 
listed in Table 1. Since the HD model (/cTin ~ 1.16 keV) is known to be inconsistent for ULXs, 
we focus on the p-free model below. 

Upon investigating the fitting carefully (see Figs. 2b and 2c), we see that the jo-free 
disk model (plus the Raymond-Smith component or not) can not provide an acceptable fit to 
the band above 7 keV. The hard tail above 7 keV shown in the plot suggests the existence 
of an additional hard component. Following previous works on BH XRBs (e.g., Kubota & 
Makishima 2004), we add to the p-free disk model an extra power-law component (Fig. 2(i), 
and the fit is significantly improved with a confidence level above 99.99999% based on the F- 
Test. Noting that the F-test was questioned in testing the significance of an additional spectral 
component (Protassov et al. 2002), we performed Monte-Carlo simulations to demonstrate the 
significance of the power-law component. In terms of the best-fitted model without a power-law 
component, we made 1000 artificial spectra, and run spectral fitting by adding an extra power- 
law component. In the 1000 simulations, the extra power-law component can only improve 
the fitting with A^^ < 16, far below the actual Ax^ = 38.65 in the real spectra, which means 
that the confidence level is far above 99.9%. Therefore, the extra power-law component is 
statistically solid, and the same statement applies to the Raymond-Smith component stated in 
the above paragraph. 

Furthermore, with the following reasons we can rule out the possibility that the hard 
tail in the spectra is due to improper background subtraction. First of all, the hard tail is not 
only visible in the co-added spectra, it is also obvious while we fit all the individual PN, MOSl 
and M0S2 spectra together. This implies that the hard tail is not due to the possible improper 
background subtraction during the spectrum co-adding. Secondly, the hard tail is also visible 
in three of the individual observations with the highest number of data bins. Spectral fitting 
to individual observations also rules out the possibility that the hard tail is dominated by one 
single observation. 

4. Discussion 

Various models have been proposed to explain the spectra of ULXs in literature. 
However, due to limited photon counts, in many cases it is not possible to distinguish these 
models by spectral fitting only. Extra constraints on these models have been given by analyzing 
the spectral evolution, as reported for several ULXs, e.g., IC 342 (Kubota et al. 2001), NGC 
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1365 X-1 (Soria et al. 2007), and Holmberg IX X-1 (La Parola et al. 2001). After investigating 
the spectral evolution of NGC 1313 X-2, Feng & Kaaret (2007) found that the variation of the 
accretion disk component deviated significantly from L cx relation if fitted with MCD+PL 
model, while roughly consistent with L oc if fitted with the p-free model. They suggested 
that this source supports the slim disk model and is against the MCD model. However, flux 
variations by more than one order of magnitude are rare in ULXs, thus such a technique can not 
be applied to most ULXs, such as M33 X-8 we studied here, which shows only small amplitude 
variations. 

In this paper we present detailed X-ray spectral fitting to 12 XMM exposures on the 
ULX source M33 X-8. We found no significant rapid variations within individual observations. 
The X-ray flux remains persistent during 10 of these 12 observations with the flux variation 
amplitude < 20%, and no signiflcant spectral evolution was detected between observations. 
Then the data of these 10 observations were combined to derive composite spectra with much 
higher S/N ratio for spectrum fltting. We flnd that the MCD-I-PL model and the p-free disk 
model provide comparable flts to both individual and composite spectra. 

However, a signiflcant hard tail above 7 keV is detected in the residual spectra of the 
p-free disk model. We added a power-law component to the p-free disk model to represent this 
hard tail, and flnd that the p-free disk -|- PL model provides the statistically best flt comparing 
to other models. The probable reason for the power-law component in this model is that the 
inner disk becomes optically thin in super-Eddington cases (Artemova et al. 2006), or there is 
a comptonized corona (Gladstone et al. 2009) similar to that of a standard thin accretion disk. 
The fluxes of the p-free disk component and the power-law component in 0.3-10.0 keV band are 
1.43x10"^^ ergs cm^^ s~^ and 0.83 xlO~^^ ergs cm"^ s~^, respectively, with a ratio of 0.63:0.37; 
this indicates that the p-free disk component dominates over the power-law component. 

The X-ray luminosity of M33 X-8 is 1.7 x 10'^^ ergs s~^ for the given distance of 0.7 
Mpc (Ho et al. 1997). The luminosity and the best-fltted inner disk temperature kTi^ = 1.38 
keV indicate that the mass of the black hole in M33 X-8 is ~ IOMq (see Fig. 1 of Watarai 
et al. 2001), in good agreement with previous works (Makishima et al. 2000; Foschini et al. 
2004). Because the X-ray emission is dominated by the disk component, the mass derived here 
should be reliable (Vierdayanti et al. 2006). We note that the luminosity exceeds marginally 
the Eddington luminosity of a ~ IOM0 black hole. In this case, the X-ray emission from the 
outer region of a standard thin disk could not be neglected, and the moderate value of p = 0.571 
is well consistent with the theoretical calculation (Watarai et al. 2000). In Galactic black hole 
X-ray binaries, the very high state (or steep power-law state) is deflned as a state in which the 
luminosity is exceedingly high (Lx > 0.2 LEdd ) and the X-ray spectrum displays substantial 
nonthermal radiation, which may constitute 40-90 % of the total flux, with a photon index 
larger than 2.4 (McClintock & Remillard 2006). The state of M33 X-8 is more like the thermal- 
dominant state of X-ray binaries in the Galaxy, but with much higher luminosity, occupying a 
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new ultraluminous accretion state (Gladstone et al. 2009). 
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Fig. 1. The 0.3 - 10.0 keV lightcurve of M33 X-8. Obs #4 and #8, which show significantly higher 
luminosity level are marked as solid circles. The solid line marks the mean luminosity of the rest 10 
observations. 120% and 80% of the mean luminosity arc marked by the dashed lines to illustrate the 
amplitude of the variation. 
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Fig. 2. Fitting to the combined spectra of ULX M33 X-8. Panels (a) and (&) show the fitting residuals 
for the HD model and the p-iree model, respectively. After adding a Raymond-Smith component to the 
models, panels (c) and (d) show the residuals for the p-bee model and the p-bcc + power-law model, 
respectively. The combined spectra and the best fitted model (the p-free -I- power-law model) are plotted 
in the upper panel. 



10 



Table 1. BEST-FIT SPECTRAL PARAMETERS OF M33 X-8 
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8 0141980301 

9 0141980401 

10 0141980501 

11 0141980601 

12 0141980801 
Combined 



PN/M1/M2 
PN 

PN/M1/M2 
PN 

PN/M1/M2 
M1/M2 
^41/1^2 
PN 

M1/M2 

PN/M1/M2 

PN 

PN/M1/M2 
PN/M1/M2 



7.4/10.5/10.5 
3.4 

1.7/9.4/9.5 
8.0 

9.0/12.1/12.1 

12.1/12.1 

8.5/9.7 

10.0 

8.3/8.6 

1.5/7.9/8.1 

10.5 

7.6/10.0/10.0 
41.0/79.0/80.8 



0-148 



1.43 
2.06 



1.50 

1.55 
1.54 
1.40 



1.56 
2.54 
1.67 



0.563 



2.26 
2.33 
2.27 
3.47 
2.23 
2.49 
2.25 
2.20 
2.70 
1.95 
1.93 
2.08 
2.18 



1240.79/1166 

88.65/84 

565.81/579 

208.85/199 

1087.61/1026 

426.00/467 

406.17/382 

238.70/241 

422.10/375 

748.04/757 

428.83/456 

1315.83/1164 

1679.66/1440 



phabs* (ray+diakhb+po) 


Combined 


PN/M1/M2 


41.0/79.0/80.8 


0. 


,00 + 0-009 
■^■^"-0.009 


161 + 0-015 1 i6 + l)-lM 
-0.016 -0.03 




2.30.0.05 


2.45 


1574.66/1435 


phabs* (ray+diskpbb) 


Combined 


PN/A'I1/M2 


41.0/79.0/80.8 


0. 


-„ + 0.U06 
-^^-^-0.005 


153 + 0020 +0.03 

o i-'^-o.oia i-^o„o 04 


535 + 0-OU'l 
0-^35_o.oo4 




2.23 


1578.18/1437 


phabs* (ray-|-diakpbb-i~po) 


Combined 


PN/M1/M2 


41.0/79.0/80.8 


0. 


i^^^^-.S?! 


lfi4+0 009 , OO + 0.09 
010*-0.013 l-3**-0.08 


„ ,7, +0.032 
O-^'^ -0.030 


2.o, + 0.11 


2.23 


1539.53/1434 



Instruments: data from which instrument. PN, MOSl(Ml) or MOS2(M2), are used; Exposure: clean exposures for 
corresponding instruments after background flares excluded; A^h- column density along the line of sight; kT: plasma 
temperature; kT-m'- inner disk temperature; p: the temperature gradient; T: power-law photon index ; fx'- 0.3-10 keV 
intrinsic flux in the units 10^^^ ergs cm~^ s^^; x^/dof: ^^^d degree of freedom for the best-fit model. 
Combined: data from combining all observations except for Obs #4 and ^8. 
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